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Abstract 

In a recent paper I have introduced a package for the exact simulation of power-law 
noises and other colored noises (E. Milotti, Comput. Phys. Commun. 175 (2006) 
212): in particular the algorithm generates l/f a noises with < a < 2. Here I 
extend the algorithm to generate l/f a noises with 2 < a < 4 (black noises). The 
method is exact in the sense that it produces a sampled process with a theoret- 
ically guaranteed range-limited power-law spectrum for any arbitrary sequence of 
sampling intervals, i.e., the sampling times may be unevenly spaced. 
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NEW VERSION PROGRAM SUMMARY 

Program Title: PLNoise 
Catalogue identifier: 

Licensing provisions: none 

Programming language: ANSI C. 

Computer: Any computer with an ANSI C compiler: the package has been tested 
with gcc version 3.2.3 on Red Hat Linux 3.2.3-52 and gcc version 4.0.0 and 4.0.1 on 
Apple Mac OS X-10.4. 

Operating system: All operating systems capable of running an ANSI C compiler. 
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RAM: the code of the test program is very compact (about 60 Kbytes), but the 
program works with list management and allocates memory dynamically; in a typ- 
ical run with average list length 2 • 10 4 , the RAM taken by the list is 200 Kbytes. 

Classification: 4.13 

External routines: The package needs external routines to generate uniform and 
exponential deviates. The implementation described here uses the random number 
generation library ranlib freely available from Netlib [1], but it has also been suc- 
cessfully tested with the random number routines in Numerical Recipes [2]. Notice 
that ranlib requires a pair of routines from the linear algebra package LINPACK, 
and that the distribution of ranlib includes the C source of these routines, in case 
LINPACK is not installed on the target machine. 

Catalogue identifier of previous version: ADXV_vl_0 

Journal reference of previous version: Comput. Phys. Commun. 175 (2006) 212 

Does the new version supersede the previous version?: Yes 

Nature of problem: Exact generation of different types of colored noise. 

Solution method: Random superposition of relaxation processes [3], possibly fol- 
lowed by an integration step to produce noise with spectral index > 2. 

Reasons for the new version: Extension to l/f a noises with spectral index 2 < 
a < 4: the new version generates both noises with spectral with spectral index 
< a < 2 and with 2 < a < 4. 

Summary of revisions: although the overall structure remains the same, one routine 
has been added and several changes have been made throughout the code to include 
the new integration step. 

Unusual features: The algorithm is theoretically guaranteed to be exact, and unlike 
all other existing generators it can generate samples with uneven spacing. 

Additional comments: The program requires an initialization step; for some pa- 
rameter sets this may become rather heavy. 

Running time: running time varies widely with different input parameters, how- 
ever in a test run like the one in section [3] in the long write-up, the generation 
routine took on average about 75 /is for each sample. 

References: 
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Recipes in C: The Art of Scientific Computing, 2nd ed. pp. 274-290 (Cambridge 
Univ. Press., Cambridge, 1992). 
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LONG WRITE-UP 



1 Introduction 

I have recently developed a power-law noise generator which uses a superposi- 
tion of uncorrelated simple relaxation processes with a uniform distribution of 
relaxation rates to produce 1/f noise, and a power- law distribution of relax- 
ation rates to obtain l/f a noise [TH2] : the generator can be easily modified to 
produce complex superpositions of different power-law noises, like those ob- 
served in precise clocks [3] . This generator also has a very unusual and unique 
feature: the superposition mechanism takes into account the correlation be- 
tween different samples, and does it so exactly, so that its noise output can 
be sampled at arbitrary sampling times. This is important whenever sampling 
cannot be performed at evenly spaced times: uneven sampling is quite common 
in astronomical observations [I], but also in other fields such as climatology 
[5], it shows up whenever observations have unintentional gaps and missing 
data [6], and sometimes bunched sampling may even be a technical need, 
e.g. in hard-disk controllers [7j; moreover, uneven sampling has some special 
properties, because it is not band limited in the same sense as ordinary, even 
sampling [8J. The generator does not produce exactly Gaussian noise, but the 
closeness to Gaussianity is tunable and the noise distribution can approach a 
true Gaussian to any required degree by a proper parameter setting [1] . 

Most of the observed power-law (1/ f a ) noises have spectral indexes < a < 2, 
with an apparent clustering around a — 1. However noises with higher spec- 
tral indexes 2 < a < 4 also show up in several unrelated systems PfTU] like the 
water level of the Nile river, economics, orchid population size [IT] and local 
temperature fluctuations and affect precise timekeeping [3] and our ability to 
predict environmental and animal population variables [12] . Noises with a > 2 
also appear in the energy level fluctuations of quantum systems [T3"1[T4"] . These 
noises correspond to generalizations of the standard one-dimensional random 
walk, and for this reason the corresponding process is often called fractional 
Brownian motion (fBm) (likewise Gaussian noises with < a < 2 are also 
called fractional Gaussian noises (fGn)). Because of their extreme peaking 
behavior at low frequencies these noises are also called "black" [10], and they 
display marked persistence properties [9] that may lead to the mistaken iden- 
tification of underlying trends in experimental data [15J. 

It is easy to see that there is no way to produce black noises from the superpo- 
sition of simple relaxation processes, the spectral density of a train of random 
pulses with simple exponential response with decay rate A is proportional to 
1/{oj 2 + A 2 ) and its derivative in a log-log plot ranges from to -2: thus no 
simple superposition can have a spectrum that decays faster than 1/f 2 and 
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the algorithm described in [T1I2] is unable to simulate this very interesting class 
of power-law noises. This paper addresses the practical problem of including 
black noises in the the noise generator introduced in [2J, preserving all its good 
features. 



2 Generation of black noises 

It may be argued that the superposition argument could still be used if one 
takes pulses with non-exponential shapes, however in that case the theory 
exposed in [1] must be revised and the analysis becomes considerably harder. 
These difficulties can be circumvented when we note that the relationship 
between the spectral density S x (u) of a random process x and the spectral 
density S y {u) of its integral y is S y (u) = S x {uj)/uj 2 and therefore it is possible 
to generate a black noise taking the integral of a power-law noise obtained from 
simple superposition. Here, as in [1112] . I consider a noise signal x(t) which is 
a linear superposition of many random pulses, i.e., pulses that are random in 
time and are associated to a memoryless process with a Poisson distribution, 
and such that their pulse response function is 



with a decay rate A which is a positive definite random variate with probability 
density g\(\) in the range {X min , X ma x) so that 



where tf. is the time at which the fc-th pulse occurs, Ak is its amplitude and 
Afc is the decay rate of its pulse response function. In principle the amplitude 
is also a random variate, but here I take a fixed amplitude = A. If n is the 
pulse rate, it can be shown pQ that the average signal level is 







k 




(3) 



the variance is 




(4) 
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and the spectral density is 
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where the decay rate A which is a positive definite random variate with prob- 
ability density g\{X) in the range (X min , X max ). I define the normalized zero- 
mean process as follows 

x(t)-(x) AJ2k h (t ~tk,h) - (x) ,„x 
x N{t) = i = , , (6) 

v / ((^) 2 ) v^) 2 ) 



and the integral of X]y(t) is 

y N (t)= [ x N (t')dt>= [ ^2=Mdt' (7) 

'((Ax) 2 ) 



. Equation (JTj) is a trivial analytical answer which is still useless for inclusion 
in the generator and it must be suitably rearranged: from equation (jTJ) we 
see that = 0, and that is a kind of continuous one-dimensional 

random walk process 



t+At 



y N (t + At)=y N (t) + J x N (t')dt' 



t 

t+At 



: y N (t) + . 1 I / x(t')dt' - (x)At I (9) 



v / ((^) 2 ) V i 



where 2/jv(t) is the integrated process sampled at time t and yn{t + At) is 
sampled at time t + At (notice that At is not a short integration time, but 
the arbitrary time interval between any two sampling times). From equation 
(J2J) we note that the integral in equation ([9]) can be rearranged as follows 



t+At 



t+At 



J x(t')dt' = AY, J h(t' - t k , X k )dt' 
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where the summations in ([TO]) run over pulses that occured before t and over 
those that occurred within the (t, t + At) time interval. Finally from equations 
(IH]) and (11 01) we obtain a formula for the integrated process y^: 



1 



-\ k At 



y N (t + At)=y N (t) + 
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(11) 



and this can be incorporated in the generator, because the summations can 
be evaluated using the events of the underlying Poisson process stored in a 
linked list as in [2] . Equation (1111) is both a practical formula for the generation 
of the integrated process y^ and a generalization of the equivalent updating 
formulas in [16] (those formulas hold only for the special case of the Ornstein- 
Uhlenbeck process and its integral, i.e. for a = 2 and a = 4): the generation of 
the integrated process can be achieved producing and maintaining a sequence 
of Poisson distributed events as in |2J and using equation ( TlTi) instead of (J2D 
to evaluate the noise process. Figures [JJ and figure [2] show a pair of xn, yN 
processes obtained with the methods described above, where the integrated 
process has a spectral index a = 3.5. Figure [3] shows the spectral density of 
the integrated process, and the numerical result is in excellent agreement with 
theory. 

Finally notice that when the input process x{t) is Gaussian (and this can be 
achieved with a proper choice of generator parameters as explained in [lj), the 
difference yN(t + At) — yN(t) is a Gaussian variate as well. 



3 Changes in the library 

The changes in the library [2] are transparent for the user, and are scattered 
in several parts of the code. As in the first version, the package contains two 
header files and two code files (plus the ranlib files that are included for 
convenience, and are redistributed according to the standard Netlib rules). 
The first header file (noise. h) contains the necessary include statements and 
the structure definitions; the second header file (noise_prototypes .h) contains 
just the prototype definitions. In the present version of the generator, the 
structure used to share information between the generation routines (defined 
in the header file noise.h) is 

struct info 
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double nt; /* transition rate */ 

double tau; /* average transition time */ 

double fillUpTime; /* fill-up time estimate */ 
double f illUpLength; /* list length estimate */ 

double average; /* signal average and standard deviation */ 
double sd; 

double meaninvlambda; /* mean value of decay time */ 



double lambdamin; /* min and max decay rates */ 
double lambdamax; 

double beta; /* beta input by the user */ 

double betaO; /* actual value of beta used in the pulse distribution */ 
double lmin; /* aux. variables */ 
double lmax; 
double dl; 
double binv; 

>; 



This header also contains the definition of N decay- #define NDECAY 20.; with 
this definition the program is quite accurate (the average relative error after 
discarding the old events is ~ 2 • 10~ 9 ), however the generation process can 
be made more accurate (and slower) with a larger value of NDECAY, or less 
accurate and faster with a smaller value. 



The first code file (list -routines . c) contains the code for the list routines: 

• Append: this is a variant of the usual Push routine; 

• Process_List: this routine processes the list to get rid of the old elements 
that can no longer influence the output signal; 

• PrintXist: this routine prints the list; 

• Response: this routine computes the noise signal for a spectral index between 
and 2; 

• IntegratedResponse: this routine computes the integrated noise signal and 
is called when the user asks for a spectral index between 2 and 4; 

• Get_List_Length: this routine returns the length of the list. 

these are internal routines, they are not meant to be called by the user, and 
are listed here for completeness; IntegratedResponse is new in this version. 

The second code file (generator. c) contains the user-callable routines. These 
routines have the same names and calling sequences as the first version of the 
generator [2J, even though they have been modified in several places. 
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4 Changes in the test program 



The user interface and the output file structure of the test program included in 
the package is the same as in the previous version [2]. The important difference 
is that now the allowed values of the spectral index a are in the range < 
a < 4. 

The output file begins with a header, which is a single line with the following 
values (separated by tabs): 

(1) dt = sampling interval; 

(2) nsamp = number of samples; 

(3) tmax = time of last sample; 

(4) noise_inf o .nt = transition rate; 

(5) noise_info.tau = average time between transitions; 

(6) noise_inf o . lambdamin = A m j n ; 

(7) noise_inf o . lambdamax = X max ] 

(8) noise_inf o .beta = /3 = a. — 1; 

(9) noise_inf o .meaninvlambda = (1/A); 

(10) noise_inf o.f illUpTime = fill-up time estimate; 

(11) noise_info.f illUpLength = list length estimate; 

(12) noise_info . average = average output amplitude; 

(13) noise_inf o . sd = standard deviation of noise signal; 

The rest of the file is a set of records, each with the following tab-separated 
values: 

(1) kk = record number; 

(2) t = actual time; 

(3) tt = time of last transition event; 

(4) listLength = list length; 

(5) signal = signal (for < a < 2) or integrated response (for 2 < a < 4); 

(6) norm_signal = normalized signal value; 

If < a < 2 signal is 

signal = ex P (M* ~ ( 12 ) 

tk<t 



and norm_signal = x^, while if 2 < a < 4 signal is the output of the 
IntegratedResponse routine 



signal 
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and norm_signal = y^. 

The code distribution also contains two Mathematica notebooks to analyze 
and display the program output; the notebook display . nb is used for spectral 
index < a < 2, while display2.nb is used for spectral index 2 < a < 4. 

The test program is used to generate the example discussed in the next section. 



5 Test run 

In the example shown below the spectral index is a = 3.5 (i.e. the program gen- 
erates l// 3 ' 5 noise, and (3 = 2.5), X min = 0.0001 and X max = 1 (the power-law 
region spans approximately 4 orders of magnitude), so that this corresponds 
very closely to the example given in [2]. And indeed, the program returns a 
signal which is the time integral of the sequence generated in the example 
given in [2J, and the underlying sequence is the same (the random number 
generator is the same as that in [2J, and uses the same seed for the random 
number sequence). 

*** PLNoise *** 

1 . terminal input of control variables 

Enter sampling time interval (dt) : 1 
Enter number of samples: 4194304 
— > Total sampling time: 4.1943e+06 
Enter transition rate: 0.1 
— > Average transition time: 10 
Enter lambda_min: 0.0001 

Enter lambda_max (0 = single relax, rate, lambdamax = lambdamin) : 1 
Enter alpha (spectral index in 1/f "alpha) : 3.5 

2. initialization 

Initialization time: 0.050000 seconds 
Noise parameters: 

— Spectral shape : 
Min decay rate: 0.0001 
Max decay rate : 1 

Beta: 2.5 (spectral index is 1+beta = 3.5) 

— Poisson process: 
Transition rate: 0.1 
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Average transition time: 10 

— Algorithmic variables: 
Fill-up time: 200000 
Fill-up length: 20000 

Mean value of decay time (<l/lambda>) : 100 

Mean list length: 200 

Signal average : 10 

Signal variance : 5 

Signal standard deviation: 2.23607 

Signal skewness: 0.298142 

Rule of thumb for Gaussianity: n<l/lambda> = 10 >= 10, noise is Gaussian 

— Internal parameters : 
lmin: 0.01 

lmax : 1 
dl: 0.99 
binv : 2 

List length after initialization: 203 

. . . initialized . . . 
. . . starting now . . . 



3. main generation loop 



XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 



Generation time: 313.780000 seconds 



4. statistics 



Statistics of generated sequence: 

Average: 10.0659 

Variance: 4.91088 

Standard deviation: 2.21605 

Skewness: 0.318761 



5 . end 



The program has been compiled and run with the same compiler and on the 
same machine as the example in [2] (i.e., the compiler gcc 4.0.1 with opti- 
mization flag -03 has been used, and the program has been run on an Apple 
Powerbook G4 - 1.5 GHz with the Mac OS X 10.4.6 UNIX flavor): a compar- 
ison between the generation times shows that in this case the integration step 
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leads to a 75% increase of generation time with respect to the unintegrated 
case. 
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200000 250000 300000 350000 400000 450000 

time (arb. units) 

Fig. 1. Normalized process x^(t) with power-law spectrum l// 1 ' 5 produced by the 
generator described in |l|2j . This record contains 2 18 = 262144 samples, and the 
generation parameters are n = 10, \ m in = 0.0001, \ ma x = 1, P = 0.5 (see [I] for a 
detailed explanation of these parameters). Time does not start from because the 
initial part of the noise record is used for initialization, and the relaxation rates A 
are given in arbitrary frequency units related to the arbitrary time units used in 
the figure. 
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Fig. 2. Plot of the process yN(t) obtained from the time integration of the process 
XAr(t) shown in figure [lj as explained in the text. The apparent global linear trend 
is characteristic of black noise, and any such trend disappears or is replaced by 
another different trend in longer records. 
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0.0001 0.001 0.01 0.1 1 
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Fig. 3. Spectral density of the process yjv(^) shown in figure [2j The linear trend must 
be removed to avoid artifacts, either by detrending or by windowing: here I used 
a Hanning window. The arrows mark the positions of the minimum and maximum 
relaxation rates \ m in and A max , the thick line is the expected average spectrum, 
corrected for the window incoherent gain, and the dashed line shows the expected 
slope of a l// 3 ' 5 spectral density. The slight upward bend at high frequency is due 
to uncorrected aliasing. 
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